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Abstract —For metrology, geodesy and gravimetry in space, 
satellite based instruments and measurement techniques are used 
and the orbits of the satellites as well as possible deviations 
between nearby ones are of central interest. The measurement of 
this deviation itself gives insight into the underlying structure of 
the spacetime geometry, which is curved and therefore described 
by the theory of general relativity (GR). In the context of GR, 
the deviation of nearby geodesics can be described by the Jacobi 
equation that is a result of linearizing the geodesic equation 
around a known reference geodesic with respect to the deviation 
vector and the relative velocity. We review the derivation of 
this Jacobi equation and restrict ourselves to the simple case of 
the spacetime outside a spherically symmetric mass distribution 
and circular reference geodesics to find solutions by projecting 
the Jacobi equation on a parallel propagated tetrad as done 
by Fuchs fT2l . Using his results, we construct solutions of the 
Jacobi equation for different physical initial scenarios inspired by 
satellite gravimetry missions and give a set of parameter together 
with their precise impact on satellite orbit deviation. We further 
consider the Newtonian analog and construct the full solution, 
that exhibits a similar structure, within this theory. 
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I. Introduction 


are slightly different and the result is the so-called 
cartwheel orbit. 

iii) A more general possibility that allows to change the 
orbital plane as well as the constants of motion is the 
helical orbit configuration. 

For the precise measurements of inter-satellite distances 
in these orbit configurations the general relativistic effects 
must be investigated and their impact on observables taken 
into account for a given accuracy level. In this work we will 
focus on such orbital configurations and refer to one satellite 
as the reference object moving on the reference geodesic. The 
orbit of the second satellite will then be modeled by means 
of the geodesic deviation equation that describes how nearby 
geodesics deviate from each other due to the geometry of 
spacetime and given initial conditions. We will first describe 
the situation in standard Newtonian gravity and then turn to 
the full theory of GR to uncover relativistic modifications. 

II. Geodesic deviation in Newtonian gravity 

In the Newtonian theory of Gravity we have as central 
equations 


When satellites follow freely falling orbits around a central 
massive object like the earth, their worldlines, i.e., their paths 
through space and time must be described by the geodesic 
equation together with given initial conditions. While for some 
(past) space missions the Newtonian theory of gravity might 
be sufficient, modem and future mission scenarios certainly 
need relativistic effects to be taken into account. Thus, the 
precise description in terms of Einstein’s theory of gravity, GR, 
becomes necessary even beyond usual Post-Newtonian (PN) 
approximations. In this work we will describe the geodesic 
deviation, with satellite missions in mind, in the full context 
of GR. 

As an example, the GFZ-NASA mission GRACE-Follow- 
On KD, (a consists of two satellites which are able to measure 
the change in their relative distance (about 100 km) with an 
accuracy of the order of 10 nm. For the orbital motion of the 
two satellites we can imagine different configurations: 


AU(f) 

d^r 


—47rGp(r), 
-VU(f), 


(la) 

(lb) 


where the first of them is the field equation that relates the 
Newtonian gravitational potential U to the mass density p 
and introduces Newton’s gravitational constant as a factor 
of proportionality. Outside a spherically symmetric object we 
obtain the gravitational potential 


U{r) = 


GM 

r 


( 2 ) 


as a solutions of ([T^. The second equation ([T§ is the 
equation of motion and describes how test particles move in 
the gravitational potential given by U. The second derivative of 
the position vector is taken w.r. to universal time t that exists in 
Newtonian gravity. We can rewrite the second equation using 
index notation and obtain 


i) Tilt the orbital plane of one satellite with respect 
to the other, but keep the constants of motion (in 
magnitude) the same. This orbital configuration is 
called a pendulum orbit. 

ii) For a second configuration, the orbital plane is 
the same, but the energy and angular momentum 


y,a ^ ( 3 ) 

The coordinates are just the usual Cartesian coordinates 
and the overdot denotes derivatives w.r. to the Newtonian 
time coordinate. The argument x in the potential denotes 
the position of the test particle. Here and in the following, 
Fatin indices a,6,... take values 1,2,3. Now, we assume a 



situation as shown in Fig. i.e., we have a reference curve 
that fulfills Thereupon, we consider a second curve 
x^{t) = X^(t) + 77 " (f) and introduce the deviation vector 77 ^. 
Thus, we have 

+ 77 ^ = -d^u{x) = -d^u{x + 77 ). (4) 

Now, we linearize the right hand side w.r. to the deviation 
U{x) = U{X + r?) = U{X) + n^daU{X) + 0{if) , (5) 

and obtain finally the deviation equation in Newtonian gravity 

= -rfd^-dbUi^X). (6) 

From (|^ we clearly see that second derivatives of the Newto- 
nian potential cause non-linear deviations. If either = 0 
(homogeneous gravitational field) or U = 0 (no gravitational 
field) the deviation equation has the solution 

77"(t)=Cit + C 2 (7) 


and the deviation vector grows only linearly in time. In 
the following, we use the Newtonian gravitational poten¬ 
tial 0 outside a spherically symmetric mass distribution 
and we change to usual spherical coordinates {x^y^z) = 
(r sin 7 ^ cos (^, r sin ^ sin ip ^ r cos t^) . We further specialize the 
reference geodesic (i?(t), ©(t), to be a circular orbit in 
the equatorial plane. Thus, we have 


R{t) = Rq = const., 

( 8 a) 

0 (t) = 00 = 7r/2, 

( 8 b) 

N Igm 
= V r 3 

( 8 c) 


where the motion is oscillating with the Keplerian period 
To = 27r/co’o = 27r^r^/(GM). For this situation the deviation 
equation ([^ reduces to three differential equations for the 
components of the deviation vector ( 77 ^^, 77 ^, 77 "^), see Q, 

d^ 77 ^ 2 7 ? 

'■oj^P = 0, 


d?rf 


df2 
d77^ 

- - ^ujorf = 0 , 

d^T?^ „ d??’’ „ 

— + 2^b,— -Q. 


(9a) 

(9b) 

(9c) 


The first equation decouples from the other two and the 
solution is given by 


rf{t) = Cs coscco^ + Cq sincco^ 5 


( 10 ) 


which is an oscillation with the Keplerian period = Tq = 
211 jw{^. In |[7l the author derived only the oscillating solutions 
for the remaining two equations. However, we extend this and 
give the general solution that yields 

rf (fy = Cl + C 2 sincco^ + C 3 cosccq^ , (Ha) 

3 

= —ccoCit + 2(C2 cosccq^ ~ C3 sincco^) + C4 . 

(lib) 

There are several possibilities to perturb the circular reference 
geodesic. One can incline the orbital plane, add a constant 
radius or cause an eccentricity in the motion. All these effects 
are due to the choice of the six parameter C^. We must have 
precisely six free parameter to set the initial position and 
velocity as starting conditions. The meaning of these parameter 
and th eir impact on the perturbed orbit will be given in section 
IV-A| in one go with those for the GR results. 



Fig. 1 . The deviation of two nearby geodesics (s) and (s) = 

77^(5). Note: the deviation vector is always defined orthogonal (as measured 
with the metric g) on the four-velocity dX^{s)/ds of the reference geodesic. 


III. Geodesic deviation equation in GR 


In the spirit of general relativity we model the spacetime 
geometry, i.e., the universe we live in (or at least the relevant 
part for our model) using a pseudo-Riemannian metric tensor 

g = g^udx^dx"" ( 12 ) 

on a four dimensional manifold M. The set {M^g} describes 
the four-dimensional spacetime and local coordinates on M 
are , /i = 0,1, 2, 3. We have an affine connection V, which 
is fully defined by the Christoffel symbols 


= (13) 

On a pseudo-Riemannian manifold we can specialize V to be 
the Levi-Civita connection and we get 

+ d^^g^x - dxgf,,.) ■ (14) 


Then, the geodesic equation that describes the motion of freely 
falling particles reads 






dx^ 
ds ds 


(15) 


and gives as solution curves the geodesics x^{s) of the 
spacetime. The affine parameter s along such a geodesic 
can be interpreted as proper time and is, thus, related 
to the reading of a clock that is transported along the 
geodesic. Now, we fix a certain (known) geodesic {X^{s)) = 
(X^( 5 ),X^(s),X^( 5 ),X^(s)) and this curve will be called 
the reference geodesic in the following. To consider a neigh¬ 
boring geodesic x^{s) in a given coordinate system we make 
the ansatz 


X^{s) = X^( 5 )+77^(5) ( 16 ) 

and define, thereupon, the deviation vector 77 ^( 5 ) that connects 
both geodesics. We assume, as sketched in Fig. the four 
velocity of the reference geodesic and the deviation vector to 
be always orthogonal to each other, 

d X^ 

(17) 

Inserting ( p^ into the geodesic equation (T5\ gives a second 
order differential equation that is quadratic in the deviation 















vector. When this deviation vector is assumed to be very 
small, we can linearize with respect to the deviation vector 
itself 77 ^( 5 ) and its derivative dr]^{s)/ds if we further assume 
small relative velocities. Thus, we obtain the so called standard 
Jacobi equation 


DVW _ dX-(5)dX-(s) 


(18) 


where the covariant derivative D/ds and the curvature tensor 
components given by 


Dr?'^(s) _ d7?^(s) 

ds ds ds ’ 


(19a) 


R^,M) = - drn^ + . (19b) 

In ( p^ we can clearly see that the curvature of spacetime in¬ 
duces a possible non-linear deviation between two neighboring 
geodesics. We have seen in the Newtonian case before that 
the same effect of non-linear deviation was caused by non¬ 
vanishing second derivatives of the Newtonian gravitational 
potential and we have, thus, an intuitive understanding of the 
role of the curvature tensor; The Riemann curvature tensor 
includes second derivatives of the metric as well. A somewhat 
detailed discussion of the equation of geodesic deviation can 
be found in standard textbooks on general relativity like m 
and It should be mentioned that if we do not linearize w.r. 
to the relative velocity dr]{s)/ds a generalized version of the 
Jacobi equation is obtained, see, e.g., la, cni, oa, na. 


These constants can be derived using the Euler Lagrange 
equations 


for the Lagrangian 


d dC 


dC 


(23) 


^ (24) 

Since the metric ( [20| ) does neither depend on the time coor¬ 
dinate t nor on the angle 0 we get for the reference geodesic 


E := A(R(s)) T(s) = const., (25a) 

L := R{s)‘^ i>(5) = const., (25b) 

where the overdot denotes the derivative with respect to the 
proper time s. The general solution of the geodesic deviation 
equation ( p^ in the spacetime ( |20| ) was given by Fuchs (91 in 
terms of first integrals of the Jacobi equation. This solution 
is, however, not applicable for circular reference geodesics 
since in this case the condition dR{s)/ds = 0 holds and 
causes singularities in the equations in this reference. Shirokov 
derived periodic solutions for the Schwarzschild spacetime and 
circular reference geodesics (H. One different possibility to 
solve the geodesic deviation equation for the case of circular 
reference geodesics in the spacetime ( [20| ) is to refer it to a 
parallel propagated tetrad and solve the resulting differential 
equation in this reference system ifT^ . We will use the results 
of this method here. 


In the following we specialize the spacetime to be spheri¬ 
cally symmetric and static, described by the metric 

g = —A{r)dt^ E B{r)dr^ + r^{d'd^ + sin^ ddd^) (20) 

and introduce the coordinates x^ =t^x^ = r, = (p. 

The angles 1 ^ and p are the usual polar and azimuthal angles 
as in spherical polar coordinates and the radial coordinate r 
is defined such that circles at a distance r have circumference 
27rr. In these coordinates the reference geodesic is given by 


The simplest model that we could use to describe the 
motion of satellites in orbit around the earth is to specialize 
to the case of Schwarzschild spacetime 

A{r) = (^1 “ > -®(^) = (26) 

and the circular reference geodesic is essentially determined 
through its radius Rq and the mass M of the earth that enters 
the metric coefficients 


X°is)=Tis),X\s)=R{s), 
X\s)=e{s),X^is) = ^is). (21) 

For a metric of the form ( [20| ) we can always, without loss 
of generality, assume that the reference geodesic is fixed in 
the equatorial plane by choosing the coordinate system to 
match this condition. This is due to the spherical symmetry 
of the spacetime exhibited in ( [^ and implies i) 0(s) = 
7r/2 = const, and ii) d&{s)/ds = 0. Since we wish to actually 
describe the motion of satellites, we must restrict to timelike 
geodesics, which describe the motion of massive test particles 
at subluminal speed. For such geodesics we can identify the 
parameter s along the reference geodesic with the proper time 
according to the normalization of the four velocity 


^ _ dX^(5)dX^(5) 


( 22 ) 


Here, we use natural units in which Newtons gravitational 
constant G and the speed of light c take the value c = G = 1. 
For such timelike geodesics in the metric we obtain 
constants of motion that correspond to the conservation of 
energy E and angular momentum L, see for example (SI. 


R{s) = Rq = const. 

= — = const. 

^{s) = i> s = ^ s =: ujs 
Rq 


T{s) = ts = 
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■ s = 


A{Rq) 1 - 2M/Rq 


(27a) 

(27b) 

(27c) 

(27d) 


We should mention that the mass of the earth in the units that 
we use is about M ^ 0.5 cm. After some lengthy calculations 
one arrives at the solution of the Jacobi equation using the 
result of Fuchs in 113 
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/(«5): 


(28a) 

(28b) 

(28c) 

(28d) 

















where the two proper time dependent function f{s) and g{s) 
are given by 


fis) = c, (^1 
+ C 4 , 




{C 2 COS ks — Cs sin ks) 


(29a) 


9{s) 


Cl [ 4 ^ 

k y k^ 


+ C 2 sin ks + Cs cos ks . 


(29b) 


The constants of motion E and L as well as the remaining 
quantities k and Rs are uniquely defined by the radius Rq of 
the circular reference geodesic 


7 ? - ^ 

-^3 - 

Uq 

2 ^ M{Ro- 6 M) 
Rl{Ro-3M) ’ 
2 ^ {R0-2MY 
Ro{Ro - 3M) ’ 
, 2 ^ MRl 
Ro-3M' 


(30a) 

(30b) 

(30c) 

(30d) 


All the other constants Ci^i = 1..6 can be used to model 
different initial conditions. 


C 4 : 


C'5,6- 


The constant C 4 simply describes an offset in the 
azimuthal and temporal components, e.g., when only 
C 4 7 ^ 0 we have a motion on the same circle as 
the reference geodesic but with a constant azimuthal 
separation 

^{s) = i.(s) + ^ C^=UJS + 

Kq 

These two constants incline the orbit w.r. to the 
reference geodesic, e.g., if only C 5 7 ^ 0 one obtains a 
circular orbit with radius Rq but with a polar angle 


0 {s) = 0 (s) + 


Cs cos Co’S 
Rq 


= — ^ 60 cos Co’S 


By choosing more than one constant unequal to zero at the 
same time, we get combinations of the described effects. We 
have named the constants Ci in the Newtonian calculation in 
the same way, such that the before mentioned effects of these 
constants are qualitatively the same. 


In the following we examine two examples for possible 
orbit deviations in some detail and give a visual representation 
in Fig. Q and Q. 


IV. Results 

A. The initial conditions 


To describe different orbital configurations we have to 
examine the meaning of the constants Ci. A proper way to 
do this is to investigate the impact of each constant separately, 
i.e., having only one of them unequal to zero. The parameter 
Ci then yield the following meaning: 


Ci: 


In the case that only Ci 7 ^ 0 we obtain again a circular 
orbit in the equatorial plane with radius 

x^{s) = r{s) = R{s) + r]^{s) 


— Rq + 


ERq ^ 2 \/R^ 


— Rq + 6 r 


We have to ensure that the perturbation 6r is small in 
comparison to the reference radius Rq by means of 
choosing Ci in a proper way. 


B. Pendulum orbits 

One possible application of the above results is to model a 
pendulum orbit, where the orbital planes of two satellites are 
inclined w.r. to each other but we keep the constants of motion 
the same. For a circular reference geodesic this means that the 
radius of the second satellite’s orbit is the same and the orbit 
is therefore circular again. To achieve this, we have chosen 
Cl = C 2 = Cs = 0. The precise choice of C 5 , Cq determines 
the line of nodes and we can include an azimuthal offset using 
C 4 to prevent both satellites from colliding. A result of this 
kind is shown in Fig. for different values of the elapsed 
proper time. The perturbed motion is again circular with 

r(s)=i?o, (31a) 

0 {s) = 7r/2 + Cs/Rq cos Co’S + Cq/Rq sin CCS , (31b) 

(/)(s)=ccs + C 4 . (31c) 

This result holds equally well in the Newtonian case, there we 
can also get a perturbed orbit like this and the result shown in 
Fig. Q is the same. 


C2,3'- 


These two constants will cause an elliptical motion in 
the r — 0 plane, e.g., if only C 3 7 ^ 0 we get a radial 
and azimuthal motion of the form 


r(s) = RqE 


ERq 




Rl 


Cs cos ks = Rq^ 6 r cos ks 


0(s) = ^(s) - ^2 - ^Cs Sinks 

= ujs — 66 sinks 


For a perturbed motion of this kind, the eccentricity e 
and the semi major axis a are then given by 

6r 

e = — , a = Rq. 

Kq 


C. Cartwheel orbits 

In this section we show the results for deviation from a 
circular reference geodesic through an eccentricity. To con¬ 
struct such solutions of the Jacobi equation ( p^ we choose 
the parameter according to 

Cl = arbitrary constant, Cs = —2Ci ^ . (32) 

rv 

This choice ensures that both objects start from the same point 
in space and the radial components of the deviation vector 
77 ^ (s) as well as its derivative dg^/ds are initially (at 5 = 0) 
equal to zero. Ci will then determine the maximal possible 
radial deviation from the circular reference orbit. Note that 
this example was considered by Fuchs in ID as well, but there 




















Fig. 2. A pendulum orbit: the deviating geodesic (black) is now inclined but circular again with the same radius as the reference geodesic (red) to keep the 
energy and angular momentum the same since these are purely determined by the radius of the circular motion. We included an azimuth angle offset via C 4 to 
prevent both objects from colliding. 


the choice of the constant Cs has the wrong sign and there are 
several misprints at the indices of the constants. The perturbed 
orbit is given by 


, , „ , Bfio Cl fiih,, 

= i?o + Sr{l — cos ks) , (33a) 

i9(s) =7r/2, (33b) 

sinks^ . (33c) 

For a motion of that kind we derive the semi major axis a that 
is simply given by 


0 ( 5 ) = ( 


■Cl 


4Rs \ 4Rs 
1 -^ 


k^ J 


k^ 


a = Rq Sr 


(34) 


and an eccentricity of 


V. Conclusion and Outlook 

Modeling satellite configurations by means of the Jacobi 
equation produces qualitatively good results in the sense that 
known effects are reproduced. However, it has to be carefully 
determined to which accuracy the results can be used compared 
with two direct solutions of the geodesic equation since the 
Jacobi equation was linearized w.r. to relative distance and 
velocity between the two objects. In an upcoming paper we 
will give results for the comparison between solutions of the 
Jacobi equation and the direct solution of the geodesics equa¬ 
tion for simple spacetimes. For higher accuracy the generalized 
Jacobi equation might be used to achieve better accuracy. If 
the distance between both satellites is small but the relative 
velocity is not, one might do the linearization only w.r. to the 
deviation itself and keep higher order terms in the derivative. 
If the relative velocity is small instead but the spatial deviation 
is not, one should do it vice versa. 


The result is shown in figure for different values of 
the elapsed proper time. We can clearly see the effect of 
perigee precession that was examined already by Fuchs for this 
specific example 03. The r and ^-motions involve different 
frequencies, i.e., u and k. For a full orbit that starts at radius 
r{s = 0) = Rq and ends at r(5 = 27v/k) = Rq the elapsed 
proper time is s = 27r/k and yields an azimuthal angle 
0(5 = 27VIk) = uj 27V/k. The difference to 27v is now called 
the perigee precession A0. Thus, we get, see e.g. O, 

A(p = - 27r = 27r - 1 ) . (36) 

fV V /v / 

Inserting uj and k that were given before in terms of the circular 
radius Rq gives the known value 


For the future import steps will be to obtain solutions 
of the Jacobi equation in more realistic, but therefore more 
complicated spacetimes as models of the real earth. The 
considered Schwarzschild spacetime is the simplest model that 
does not include the rotation or even higher multipole moments 
of the earth. For these more realistic spacetimes the solution of 
the Jacobi equation will certainly need numerical integration 
or approximation methods. 
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However, in the Newtonian case there is no such quantity as 
k and only the Keplerian frequency ccq appears instead of k. 
Thus, for the Newtonian case we get no perigee precession 
and the behavior is qualitatively different in comparison to 
that shown in Fig. © - the perturbed orbit, i.e., the ellipse 
just remains unchanged. 
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